/**
 * This file is part of gapro library.
 * 
 * gapro library is free software: you can redistribute it and/or modify it under the terms of the GNU Lesser General
 * Public License as published by the Free Software Foundation, either version 3 of the License, or (at your option) any
 * later version.
 * 
 * This program is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY; without even the implied
 * warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more
 * details.
 * 
 * You should have received a copy of the GNU Lesser General Public License along with this program. If not, see
 * <http://www.gnu.org/licenses/>
 * 
 * Copyright 2012 Vinay Penmatsa. All rights reserved.
 * 
 */
package org.stupidwhiz.gapro.core;

import java.util.Arrays;

/**
 * Matrix operations with double types
 * 
 * @author Vinay Penmatsa
 * @since 1.0
 */
public class Matrix {

    double[][] matrix;

    /**
     * <pre>
     * Creates a matrix from the double type 2d array
     * Each 1d array should be a column vector.
     * 
     * So, given {{1,2,3}, {4,5,6}}
     * the matrix internally would be 
     * 
     *  --      --
     *  |  1  4  |
     *  |  2  5  |
     *  |  3  6  |
     *  --      --
     * </pre>
     * 
     * @param matrix
     */
    public Matrix(double[][] matrix) {
        this.matrix = matrix;
    }

    /**
     * Returns the number of rows
     * 
     * @return row count
     */
    public int rowCount() {
        return matrix[0].length;
    }

    /**
     * Return the number of columns
     * 
     * @return column count
     */
    public int columnCount() {
        return matrix.length;
    }

    /**
     * Returns a new matrix that is the transpose of this matrix.
     * 
     * @return the tranpose of this matrix
     */
    public Matrix transpose() {
        int m = columnCount();
        int n = rowCount();
        double[][] transpose = new double[n][m];

        for(int i = 0; i < n; i++) {
            for(int j = 0; j < m; j++) {
                transpose[i][j] = matrix[j][i];
            }
        }

        return new Matrix(transpose);
    }

    /**
     * Returns a new matrix that is the inverse of this matrix
     * 
     * @return the inverse of the matrix
     */
    public Matrix inverse() throws Exception {
        int m = columnCount();
        int n = rowCount();

        if(m != n) {
            throw new Exception("Matrix is not square");
        }

        Matrix im = identity(m);
        Matrix augmented = augment(this, im);
        reduce(augmented);

        return augmented;
    }

    private void reduce(Matrix augmented) {
        int row = 0;
        for(int col = 0; col < augmented.columnCount(); col++) {
            int k = findMaxInRowFromIndex(augmented.matrix[row], row);
            if(k == -1)
                continue;
            if(k != row) {
                exchangeRows(augmented, k, row);
            }
            double m = augmented.matrix[col][row];
            multiplyRow(augmented, row, m);
            reduceRow(augmented, row, col);
            row++;
            if(row == augmented.rowCount())
                break;
        }
    }

    private void multiplyRow(Matrix mx, int row, double m) {
        for(int i = 0; i < mx.columnCount(); i++) {
            mx.matrix[i][row] /= m;
        }
    }

    private void reduceRow(Matrix mx, int row, int col) {
        for(int r = 0; r < mx.rowCount(); r++) {
            if(r == row)
                continue;
            double m = mx.matrix[col][r];
            for(int c = 0; c < mx.columnCount(); c++) {
                double mv = m;
                mv *= -1;
                mv *= mx.matrix[c][row];
                mx.matrix[c][r] += mv;
            }
        }
    }

    private int findMaxInRowFromIndex(double[] row, int index) {
        double max = row[index];
        int k = index;
        for(int i = index; i < row.length; i++) {
            if(Math.max(row[i], max) != max) {
                max = row[i];
                k = i;
            }
        }
        if(max == 0)
            k = -1;
        return k;
    }

    private void exchangeRows(Matrix mx, int k, int row) {
        for(int j = 0; j < mx.columnCount(); j++) {
            double temp = mx.matrix[j][k];
            mx.matrix[j][k] = mx.matrix[j][row];
            mx.matrix[j][row] = temp;
        }
    }

    /**
     * Creates a n x n identity matrix. Elements are initialized to 0 and diagonal elements to 1
     * 
     * @param n
     * @return identity matrix
     */
    public static Matrix identity(int n) {
        double[][] elems = new double[n][n];

        int key = 0;
        for(int d = 0; d < elems.length; d++) {
            Arrays.fill(elems[d], 0);
            elems[d][key] = 1;
            key++;
        }
        return new Matrix(elems);
    }

    public static Matrix augment(Matrix m, Matrix n) throws Exception {
        int ml = m.rowCount();
        int nl = n.rowCount();

        if(ml != nl) {
            throw new Exception("Row count of m, n must be equal");
        }

        int colCount = m.columnCount() + n.columnCount();
        double[][] augmented = new double[colCount][ml];

        for(int i = 0; i < m.columnCount(); i++) {
            System.arraycopy(m.matrix[i], 0, augmented[i], 0, ml);
        }

        for(int i = 0; i < n.columnCount(); i++) {
            System.arraycopy(n.matrix[i], 0, augmented[i + m.columnCount()], 0, nl);
        }

        return new Matrix(augmented);
    }

    /**
     * Returns the string representation of the matrix as column vectors
     * 
     */
    public String toString() {
        StringBuilder sb = new StringBuilder();
        for(int j = 0, m = columnCount(); j < m; j++) {
            sb.append("<");
            for(int i = 0; i < rowCount(); i++) {
                sb.append(matrix[j][i] + (i == rowCount() - 1 ? "" : ","));
            }
            sb.append(">");
            if(j != m - 1)
                sb.append(" ");
        }
        return sb.toString();
    }
}
